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Abstract 

Using the complex Langevin sampling strategy, field theoretic simulations are performed to 
study the equilibrium phase behavior and structure of symmetric polycation-polyanion mixtures 
without salt in good solvents. Static structure factors for the segment density and charge density are 
calculated and used to study the role of fluctuations in the electrostatic and chemical potential fields 
beyond the random phase approximation. We specifically focus on the role of charge density and 
molecular weight on the structure and complexation behavior of polycation-polyanion solutions. A 
demixing phase transition to form a "complex coacervate" is observed in strongly charged systems, 
and the corresponding spinodal and binodal boundaries of the phase diagram are investigated. 
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I. INTRODUCTION 

Since statistical field theory was first applied to self avoiding polymers by Edwards 
polymer field theory models and techniques have been developed and refined for a wide va- 
riety of polymeric systems [3, 3, Q, E, Q, zj|- In the field theoretic approach, one often has to 
invoke approximations in order to make functional integrals for partition functions and av- 
erage properties tractable. A convenient startingpoint is the mean field approximation, also 
known as self- consistent field theory (SCFT) |l|Jj]. Since the effective coordination number 
grows as the square root of the molecular weight in concentrated polymeric systems, SCFT 
is argued to be asymptotically exact in polymeric melts when the degree of polymeriza- 
tion becomes infinite. Moreover, the theory serves as the reference for more sophisticated 
approximation schemes that attempt to account for field fluctuation effects 2|, |8|. Over 
the last decade, with the ever increasing power of digital computation, numerical SCFT 
has become a routine tool for investigations of the structural and thermodynamic proper- 
ties of inhomogeneous polymers, including polymer alloys and block copolymers of varying 



architecture 



10 



11]. 



Despite the success and the popularity of SCFT among polymer researchers, there are 
classes of systems where SCFT is known to be inaccurate j8(. These systems, character- 
ized by strong or non-negligible density fluctuations, include polymer solutions in dilute 
and semi-dilute regimes, systems close to a phase transition, and polyelectrolyte solutions, 
among others. The mean field approximation (i.e. SCFT) is particularly ill-suited to poly- 
electrolytes; indeed, for some charged polymer systems the Coulomb interaction does not 
contribute to the free energy at the SCFT level. Such a system is a homogeneous poly- 
electrolyte solution of arbitrary concentration in the bulk. The overall charge neutrality 
condition of the system leads to a constant self-consistent electrostatic potential field that 
makes no contribution to the free energy of the solution. Hence, in such a system the electro- 
statics contributes to the free energy only via positional correlations among charges. Such 
charge correlations can be very strong in polyelectrolytes due to the high valency and the 
low translational entropy of macroions. 

This strong charge correlation effect in polyelectrolytes was noticed nearly a century 
ago {12I ]: mixtures of oppositely charged biopolymers or synthetic polyelectrolytes can yield 
dense liquid precipitates coexisting with supernatant solvent under standard physiological 
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conditions at room temperature. This liquid-liquid phase separation, referred to as complex 
coacervation, is a manifestation of charge correlations and fundamentally differs from the 
macrophase separation that typically occurs in solutions or melts of incompatible neutral 
polymers. Applications of complex coacervates exist in both nature and technology. As an 
elegant example of the former, the "sand castle" worm constructs habitats on the meter 
length scale by gluing millimeter-sized grains of sand together under the sea using a wonder 
glue composed of anionic and cationic proteins jl3| . Technological applications span the 
fields of water purification, adhesives, coatings, and biotechnology. For example, DNA 
sensors are being developed based on the complexation of target DNA molecules (anionic) 
with synthetic conjugated cationic polyelectrolytes 14]. Other potential applications in 
biotechnology relate to drug delivery and gene therapy, and invoke charge complexation to 
build useful structures with a payload (such as DNA) for delivery in an aqueous environment. 
An example of such a structure is the self-assembled polymeric micelles produced by mixing 
solutions of polyanions and polycations, where one or both polyelectrolytes contains a charge 



neutral hydrophilic block 15, 



16 



17 



18j |. The core of each micelle is a coacervate composed 



of negatively and positively charged polymer segments, while the corona consists of water 
soluble neutral blocks that can serve to protect the payload in the core. 

Theoretically, the thermodynamics of complex coacervation was first examined by Over- 



beek and Voorn 



m 



20J. Although they correctly identified the driving force for the phase 



separation as a Coulomb attraction between oppositely charged macroions, they neglected 
the charge connectivity along the polymer backbone and relied on the simple Debye-Hiickel 
approximation in calculating the electrostatic free energy, which is valid only in very dilute 
systems [20]. Since that time, the main issue in the theory of complex coacervation has been 
how to properly calculate the electrostatic charge correlation in polyelectrolyte solutions. 
Following the maturation of polymer field theory, significant progress in our understand- 
ing of charge correlation phenomena have resulted from application of the random phase 



approximation (RPA) to weakly charged polyelectrolytes 21 



22 



23 



26, 



24J. and loop expan- 



271 ] . Unfortunately 



sions to treat stronger correlations beyond the level of the RPA [25l . 
these analytical techniques are still limited in scope and are particularly difficult to apply 
to inhomogeneous structures, such as mesophases. 

The pervasiveness and the importance of polyelectrolyte systems in nature and biological 
applications, coupled with the obvious inadequacy of SCFT in describing those systems, 
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motivates the development of a systematic numerical way of incorporating field fluctuation 
effects beyond SCFT. Ideally, such a scheme should be capable of treating an arbitrary field 
theory model in the absence of any approximations (aside from numerical errors that arise 
from resolving and statistically sampling the fields), and thus be capable of describing strong 
charge and density correlation effects. Beyond serving as a general purpose simulation tool, 
such a "field theoretic simulation" method could be used as a test bed to validate analytical 
results based on loop expansions and other approximations. 



Our group has developed a broad suite o 
of polymer statistical field theory models 



methods for conducting numerical simulations 



281 ] . Since the models typically have effective 



Hamiltonians that are complex, rather than real, phase oscillations thwart conventional nu- 
merical simulations based on Monte Carlo sampling. We have found that this "sign problem" 
can be effectively treated by adopting the complex Langevin (CL) sampling method from 



nuclear physics 



29 



30] , and applying it in tandem with advanced numerical methods. This 



emerging field of field theoretic simulations (FTS) of polymeric systems has already been 
applied to several polymer models where SCFT is known to be inaccurate, such as block 
copolymers near their order- disorder transition (ODT) 3lU32|. 



shitz point 



ternary blends close to a Lif- 



33] , and polymer solutions in the semidilute regime 



34 



351 ] . The FTS technique, 



however, is generic enough to be applicable to a much broader class of field theory models 



28 



36]. 



of complex fluids, and is potentially expandable to systems out of equilibrium 

In this article, we report on the application of FTS to polyelectrolyte systems. Based on a 
field theoretic model of a simple binary polyelectrolyte mixture, FTS is used to study density 
correlations in the system and to directly monitor the complexation phenomena that leads to 
the formation of a coacervate. We recently contributed a short highlight article that included 
some preliminary findings [37]. Here we provide a more complete set of results, a detailed 
discussion of our model and numerical methods, and an in-depth analysis. A particular 
focus of the present paper is on charge and density correlations in the solution. Structure 
factors of segment density and charge density are calculated using FTS and compared with 
RPA structure factors. 
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II. THEORY AND NUMERICAL METHODS 



In this section, the field theory of our model polyelectrolyte solution is introduced. The 
thermodynamics of the system in the mean field approximation is studied, followed by a 
discussion of the first non-vanishing (RPA) correction to SCFT assuming small amplitude 
field fluctuations. Next, the FTS technique is described along with CL sampling. The 
section ends with a discussion of the numerical methods used to integrate the CL equations. 

A. Field theory model 

Here we propose a simple, yet fundamental, model that is capable of giving rise to the 
phenomenon of complex coacervation with a minimal set of parameters. Specifically, we 
consider a solution consisting of a mixture of polycations and polyanions in an implicit good 
solvent with a uniform dielectric constant e. For simplicity, we choose to work with a sym- 
metric model in which the two types of polyelectrolytes are identical except for the sign of 
the charge that they carry. Half of chain molecules in the system are positively charged 
and the other half are negatively charged, thereby preserving the condition of electroneu- 
trality. Again, for simplicity, we do not include either counterions or salt, although such 
generalizations are straightforward. Our highly idealized model system, however, could be 
approximately realized by mixing a polyacid with a polybase of equal molecular weight and 
equal but opposite charge in water. We formulate our model in the canonical ensemble and 
consider a mixture of n polyanions and n polycations in a volume V of a three-dimensional 
physical space. The polymer backbones are modeled as discrete Gaussian (bead-spring) 
chains with degree of polymerization N; the number of beads per chain is N + 1. In addi- 
tion, we also employ the continuous Gaussian chain model to derive some analytical results 
in section III C( this model is a large- N limit of the discrete Gaussian chain. The charge per 
bead is +z for the cation species and —z for the anion species; all other physical charac- 
teristics of the two types of chains are assumed to be identical. Interactions among beads 
include intramolecular spring forces, and intra- and inter-molecular excluded volume and 
Coulomb interactions. 

The classical canonical partition function of the model just described can be expressed 
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(in "particle" form) as 



2rc 



Z c (2n, V, T) = 2 \ 2n(JV+1) ( J] I dR a I exp[-/3C/ - (3U ev - es ] , (1) 

(n!) (A T ) \q=i^ / 

where At is the thermal de Broglie wave length for unconnected beads and (3 = l/(fcsT) is 
the reciprocal of the thermal energy. In indexing chains with a, the first n chains are the 
polyanions (1 < a < n) and the remaining n chains are the polycations (n + 1 < a < 2n). 
The chain conformation vector R a is a set of coordinates of beads belonging to chain a: 
R-a = {i"a,o, r a,i, *a,2, • • ' , Toi,n}- The set of bead coordinates for all 2n chains is denoted by 
R. j3Uo is the intramolecular, short-ranged potential that connects beads along the chain 
backbone: 

2n N 

(3Uo[K] K\v aJ - raj-il). (2) 

«=1 j=l 

Gaussian discrete chains are assumed with harmonic spring potential (3h(r) = 3r 2 /26 2 , where 
b is the length of each statistical segment. 

The remaining contributions to the potential energy are potentials of mean force (denoted 
by an over bar) because the solvent is treated implicitly. (3U ev is the excluded volume 
interaction of the system: 

ev [R] = ^Jj dhdh' pn(T)u(\T - r'\)p n (r'). (3) 

Here the microscopic segment (bead) number density operator p„(r) is defined as the sum 
of the microscopic polyanion segment density p_ (r) and the microscopic polycation segment 
density p + (r), where 

n N 

P-(r) = J2J2 S ^- r ^) (4) 

o=l j=0 

and 

2n N 
a=n+l j=0 

We also adopt Edwards' simple delta function model for the volume interaction [l|: 

Pu{\r\)=u 6(r), (6) 
where uq is the excluded volume parameter. 



6 



The final contribution to the energy, /3U es , is the electrostatic interaction between charges 
in the system: 

PU es [R] = l -JJ ^Vp c (r)^-^p c (r') , (7) 

where Ib = /3e 2 /e is a constant Bjerrum length; e is the charge of a proton. Since every 
bead carries charge +z or —z, the microscopic charge density is defined as 

Pc(r) = -zp-(r) + zp + (r) . (8) 

By electroneutrality, the volume integral of this microscopic charge density vanishes. 

Our next step is to utilize Hubbard- Stratonovich transformation to convert the "particle" 
representation of the partition function in Eq. (CEJ) to a more convenient statistical field theory. 



In this process, two auxi 



be recast in the form 



28] 



iary fields, and w, are introduced and the partition function can 



Z C = Z jvw J X?0 exp (-H[w, 0]) . (9) 

where the functional integrals over the two fields are taken in the real function space and 
H[w, 0] is a complex effective Hamiltonian. The field w(r) can be interpreted as a fluctuating 
chemical potential field conjugate to p n (r), while 0(r) is an electrostatic potential field con- 
jugate to p c (r)- Zq is a prefactor that includes the ideal gas entropy of non-interacting chains 
and some spurious self-interactions contained in Eqs. ([3]) and (171). These self- interactions 
produce only a constant shift in chemical potential and have no thermodynamic consequence. 
The effective Hamiltonian corresponds to the functional 

H[w,<j>] = - [ d 3 r ( l V< ft( r )l | - n ]nQ\iw+iz(f)]-n\nQ\iw-iz<j)], (10) 

2 J \ u AtxIb J 

where Q[iw ± iz(j)\ is a single chain partition functional of decoupled poly electrolytes in the 
conjugate fields. Specifically, Q is defined as the ratio of the partition function of a single 
chain subject to the (pure imaginary) fields iw(r) ± iz<j)(r) to the partition function of an 
ideal chain, 

Of +',A1 Jd 3 r N+1 eM-(3Us(r N+1 )] nn 
Q hw ± tzo \ — 7T , (11) 

v(fd*heM-PH\m 
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where r N+l denotes the bead coordinates of the chain and 

Af N 

(3U s (r N+1 ) = - r^l) +i^[ W (r i )±#(r i )]. (12) 

The functional Q[iw ± izcj)] is thus normalized so that Q[0] = 1. In practice, Q is com- 
puted for an arbitrary field ip(r) according to Q[ip) = y J dr 3 q(r,N; [ip]), where q(r,j; [ip]) 
represents the statistical weight for the jth bead of a chain to be at position r. This ob- 
ject g(r, j; [■?/>]), referred to as a chain propagator, is calculated by iterating the following 
Chapman-Kolmogorov type equation 

q(r,j + lM) = ' exp[-^(r)] J dh' q{v' , j- [V>]) exp (- ^ (13) 

from the initial condition of q(r, 0; [ip]) = exp[— ip(r)]. It is notable that the integral on 
the right hand side of Eq. (fl3j) is of convolution form so can be efficiently evaluated using 
Fourier transforms. 

The average of any thermodynamic observable Q can be formally defined as an ensemble 
average of the corresponding operator Q[w,(j)] over the auxiliary field variables 

In the present paper, we are particularly interested in operators for densities of polymer 
segments and charges, and operators whose averages over the field variables yield two- 
point density and charge correlation functions. The latter can be obtained by augmenting 
Eq. with a source term involving external fields conjugate to the microscopic segment 
and charge densities, and using In Zc as a generating functional for the cumulant moments 



of density 



28j |. The segment density operator for polycations (polyanions) is found to be 



b\wQ\iw ±iz(t)} 

P±{r;[w,(p\) = -n , , 15 

o(tw(r) ± tz<p(r)) 



so that (p±(r)) = (p±(r; [w,<p])). In practice, the functional derivative on the right hand 
side of Eq. ( T15I) can be computed using the chain propagator q as 

The total segment number density operator p n (r) is the sum of p+(r) and p_(r): 

Pn(r; [w, 4>}) = p+(r; [w, (j)}) + p_(r; [w, <j>]). (17) 



Likewise, the charge density operator p c (r) is the sum of p+(r) and p_(r) weighted by the 
respective charge densities: 

Pc(r; [w, 0]) = z {p+(r; [w, <p\) - p_(r; [w, 0])} . (18) 

The correlations in the density fluctuations can be formally related to the correlation 
functions of the corresponding auxiliary fields. The pair correlation function of total segment 
number density can be computed from the pair correlation function of the auxiliary chemical 
potential field according to 

(A.(r) W 0> = ^^-^^>. (19) 

The charge density correlation function can be similarly related to the pair correlation 
function of the auxiliary electrostatic potential field 

The cross correlation between the segment number density and the charge density is similarly 
given by 

Mr)V 2 0(r')) 

{pn (r )p c r = — . 21 

UqAtiIb 

The microstructure of the polyelectrolyte solution can be characterized with static struc- 
ture factors, which are related to the pair correlation functions by Fourier transforms. The 
number density structure factor, defined as 

Snn(k) = ^jj ex P H k ■ " r')Kpn(r)pn(r')> , (22) 

can be evaluated formally from Eq. (1191) in terms of the w field pair correlation function 

1_ _ (w(k)w(- k)) 
u u 2 V 

for k ^ 0, where w(k) is the Fourier transform of w(r). Likewise, the static structure factor 
for the charge density can be calculated from the field pair correlation function as 

k* ( k> \ 2 (0(k)0(-k)) 

where 0(k) is the Fourier transform of 0(r). The structure factor for the cross correlation 
between the segment number density and the charge density can be calculated from the 
correlation between the w(r) field fluctuation and the 0(r) field fluctuation, 



Snn(k) = -- ' WW 7J " ( 2 3) 
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B. Self consistent field theory 



Self-consistent field theory (SCFT) is derived by assuming that the functional integrals 
of Eqs. ([H]) and (|14p are dominated by "mean field" configurations w*(r) and 0*(r) that cor- 
respond to saddle points of the effective Hamiltonian. The saddle point conditions produce 
the following SCFT equations: 



6H 

5w 



0=0* 



and 

SH 



— ^ + *p n (r; [w*, 0*]) = (26) 
u 



V 2 0*(r) 

— 4 J 1 1 + Vc(r; K, 0*]) = 0. 1 27) 



0=0* 



The electrostatic mean field Eq. ( 1271) recovers the conventional Poisson equation with i<p* 
interpreted as the electrostatic potential. Indeed, the physically relevant solutions of these 
equations correspond to w* and 0* being pure imaginary fields. In an unbounded system, 
or a system with periodic boundary conditions imposed on both w(r) and 0(r), and under 
good solvent conditions uq > 0, these equations yield only homogeneous saddle fields of 
constant (imaginary) w* and 0*. This trivial mean field solution corresponds to V0* = 
and w* = —iuopo, where po is the average segment number density, 2n(N + 1)/V. 

At the SCFT level of description, our model system is thus a structureless polymer so- 
lution with constant segment density and a trivial structure factor, S nn (k) = u^ 1 . The 
overall electro-neutrality makes the charge density vanish locally, because the constant pos- 
itive charge density field is compensated by the constant negative charge density field. The 
Helmholtz free energy A thus involves only the ideal gas translational entropy and the ex- 
cluded volume interaction, and is minimal when polymer chains are evenly distributed in 
the system: 

/3A = -In Z c w -\nZ + H[w*,(/)*] 

= (3A + | Po V , (28) 

where /3Aq = — In Zq is the ideal chain (translational entropy) term. Therefore, there is no 
contribution of the Coulomb interaction to the mean field free energy, and SCFT is unable 
to predict the formation of a complex coacervate. 
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C. Gaussian fluctuations 



The phenomenon of coacervation can be traced to the presence of charge correlations in 
polyelectrolyte mixtures - correlations that are neglected in the mean field approximation. 
This correlation effect can be treated analytically with the systematic loop expansion scheme. 
The first correction to SCFT in such a loop expansion can be calculated analytically for our 
model system under the assumptions of continuous (rather than discrete) Gaussian chains 
and weak (low amplitude) field fluctuations. Later, the correlation functions calculated in 
this subsection will be compared with FTS results for discrete Gaussian chains with varying 
degrees of polymerization N and for arbitrary strengths of field fluctuations. 

For the case of a polymer that is experiencing an arbitrary field ^>(r) that fluctuates only 
weakly from its homogeneous mean field ifj*, the single chain partition function Q[ip] for a 
continuous Gaussian chain can be approximated up to the second order in field fluctuation 



as 



28] 



Q[ip] «exp[-JV^(0)/V] 

N 2 
2V 2 



x { 1 + £72 E 9D{k 2 R 2 g )m^) + 0{ft) \ , (29) 



where tp(k) is the Fourier transform of VK r )> an d R g = is the radius of gyration 

of a Gaussian chain without interactions. The Debye function gr>{k 2 R 2 ) is the scattering 



function of an ideal continuous Gaussian chain 



38): 



2 

g D (x) = — [exp (-x) + x - 1] . (30) 



x 2 

With Eq. (j29p . the effective Hamiltonian of Eq. (jlOD can be approximated as 

+ 2V E [%uk)w{^(-^) + , (31) 

k^O 

where fluctuations of w and fields are decoupled at quadratic order with expansion coef- 
ficients of 

% w (k) = - + PoNg D (k 2 R 2 g ) (32) 

and 

%m = ^f B + PoNz 2 g D (k 2 R 2 g ). (33) 
11 



Using Eqs. dHJ) and (13TI) . the osmotic pressure p of the solution with weak Gaussian field 
fluctuations can be approximated as 

fd(3A\ fdlnZ c 



(3p 



Po UqPo 
N 2 



n,f3 
2 1 



dV 



24tt 



n,0 
3/2 



( 12u oPo y z 1 ( 48nl B z 2 p \ 
\ b 2 J V2 \ b 2 J 



3/4" 



(34) 



where the last negative term with square brackets is the correction to the mean field osmotic 
pressure due to field fluctuation effects. Eqn. (134]) can be expressed in a dimensionless form, 



PpR g 



C 



BC 2 



1 



2 24tt 

with reduced (dimensionless) variables of 

u N 2 



(2BC) 3/2 + 



(2EC) 



3/4 



V2 



(35) 



C 



2nRl 
V ' 



Rg 3 



E 



Aixl B z 2 N 2 
R„ 



(36) 



where C is a reduced chain concentration, B is a reduced excluded volume parameter, and 
E is a reduced Bjerrum length. 

The comparison between contributions from the mean field and the correction due to 
field fluctuations in Eq. fl35|) provides a criterion that can be used to assess the importance 
of field fluctuations. For the excluded volume interaction, the w field fluctuation can be 
regarded as weak correction when 



BC 2 > (BC) 



3/2 



(37) 



Therefore, in a three dimensional physical space, the mean field result becomes asymptoti- 
cally exact when C/B ^> 1 . In contrast, the electrostatic <fr field fluctuation contribution, 
namely the term scaling like — (EC) 3 ^, is of paramount importance because there is no 
contribution from the mean electrostatic field to the osmotic pressure. Thus, electrostatic 



1 While the condition C/B 3> 1 is commonly identified as the definition of the dense regime for neutral 
polymer solutions [28|,l3J|, we note that the precise account of the numerical prefactors in Eq. (p?I)j) produces 
a slightly different condition: C/B ^> l/18n 2 « 1/178. The numerical factor 1/187T 2 lowers the formal 
boundary between the dense and the semi-dilute regimes by more than two orders of magnitude. For 
this reason we call the case of C = 1 and B = 12 (discussed below) as "moderately dense" rather than 
"semi-dilute" , as this case may be interpreted as being on either side of the formal boundary between the 
dense and the semi-dilute regimes depending on the exact condition applied. 
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correlation effects can be neglected only when (EC) 3 / 4 is small compared with all other con- 
tributions to the osmotic pressure. Because the electrostatic mean-field term is identically 
zero, a criterion similar to Eq. (1371) cannot be derived from Eq. (135]) . and hence the range of 
validity of the electrostatic term is unknown at the one-loop (Gaussian) level. Evaluation 
of higher order terms in a loop expansion would be required to clarify the range of validity 
of this term. 

It is important to note that Eq. (1351) predicts a negative (attractive) contribution of 
electrostatic correlations to the osmotic pressure. T 



21 



22 



lis contribution, which is similar to 



can drive com- 



23 



24 



expressions derived previously using the RPA 
plexation of polyanions and polycations to produce a complex coacervate phase. As will be 
discussed below, the polymer concentration in the coacervate can be estimated by balancing 
the repulsive excluded volume terms and the attractive electrostatic correlation terms in 
Eq. (EH}. 

Using auxiliary field pair correlations calculated with the quadratic Hamiltonian of 
Eq. (13TT) . structure factors of the segment density S nn (k) and charge density S cc (k) can 
be approximated by expressions involving the quadratic expansion coefficients of Eq. fl32|) 
and Eq. (|33l) fl: 

Snn(k) « - - 2 } ... (38) 

and 2 

Scc(k) "4" %Jk)- m 

These structure factors are commonly referred to as RPA structure factors and are applicable 

in the weak field fluctuation limit of Eq. (131!) where the harmonic fluctuations in the w and 

(j) fields are decoupled. Indeed, the structure factor for the cross correlation between the 

number density and the charge density S nc (k) vanishes in this level of description, as the 

approximated Hamiltonian of Eq. fl3"TT) does not involve a term proportional to w(k)<f)(— k). 

Complex coacervation is a phase transition that occurs when the charge density correlations 

are strong enough to influence the segment density distribution. Therefore, it would seem 

that structure factors must be calculated beyond the RPA level to accurately describe the 

complexation process [25]. 
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D. Field theoretic simulations 



The analytic treatment of field fluctuation effects in the previous subsection is based on 
the weak inhomogeneity approximation of Eqs. f|2U|) and (T3TTT) . Here we turn to a direct 
numerical approach (FTS) that avoids any assumption of weak field fluctuations. 

The essence of FTS is to devise an efficient numerical strategy of multi-dimensional 
integration whereby thermodynamic averages defined in Eq. ffl4"j) can be evaluated. The 
conventional Monte Carlo importance sampling strategy is problematic for that purpose. 
Although the functional integrals in Eq. ffl~4j) are over strictly real functions w(r) and 0(r), 
the effective Hamiltonian H [w, 0] is a complex functional. Therefore, the probability weight 
proportional to exp (—H[w,<j)]) is not positive definite unless the sampling trajectory hap- 
pens to be along a constant phase path of exp(—H[w,<f>]). Because the identification of 
constant phase paths in high dimensional function spaces is computationally unfeasible, a 
fundamentally different simulation technique is required. 

The complex Langevin (CL) sampling strategy addresses the so-called "sign problem" 
associated with the non-positive definite statistical weight of the field theory 29|, |3Cj . The 
idea behind this method is to extend the real fields into the complex plane and to compute 
ensemble averages of observable quantities by sampling fields along a stationary stochastic 
trajectory in the complex function space. By extending the real fields w(r) and 0(r) into 
the complex functions W(r) = w(r) +iwi(r) and $(r) = 0(r) + z0/(r), the ensemble average 
of Eq. ( fl4|) can be reexpressed as 



(G)= / / / / VwVw/DQVfa Q[W, $}P{w, wj, 0, <j>j\, (40) 




where P[w, wj, 0, 0/] ls a rea ^ non-negative statistical distribution of fields replacing the com 
plex weight of Eq 
degrees of freedom 



4j) . This comes at the cost of doubling the number of the configurational 

m 




VwVwiV(pV(j)i S[w' — w — iwi]S[<f>' — — i(f>i]P[w, wi, 0, 0/], (41) 



Z c 

where w' and 0' are real fields. Although necessary and sufficient conditions for the existence 



of P satisfying Eq. (l4"Ti) given a complex effective Hamiltonian H have been identified 39|, |40] , 



these conditions are difficult to verify in the highly nonlinear and nonlocal Hamiltonians of 



polymer field theory. However, the current and the previous success 



CL method to polymer models provide such a proof a posteriori 
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34j, 



ul applications of the 
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4l|. 



The CL dynamics is a stochastic Langevin dynamics in the complex function space de- 
signed to generate a stationary Markov sequence of complex functions with distribution 

P[w,u> 7 ,0,0/]: 



d_ 

Of 
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and 
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4tt/ 



+ zp c (r,t; [W,$]) 



+ ^(r,i), 



(43) 



where r/ w (r, £) and ^(r, £) are rea/ Gaussian white noise fields with zero mean and variances 
proportional to the real dissipative coefficients \ w and A^, respectively, consistent with 



the familiar fluctuation-dissipation theorem of Brownian dynamics 28]. The CL equations 
should be interpreted as a fictitious, rather than physical, dynamics to sample the field 
configuration space. In the absence of the random forces, the above equations reduce to 
deterministic equations that have saddle point field configurations of Eq. (1261) and Eq. (I2T1) as 
steady state solutions. With the forcing terms, the stochastic dynamics drive trajectories in 
the complex function space that produce steady distributions P[w, Wi, (ft, <fii) that are peaked 
at saddle point configurations. Because the random forces are strictly real, the imaginary 
components of the CL dynamics drive the field trajectories towards constant phase paths, 
while the real components of the equations are responsible for stochastic motion along a 
path. Under conditions where a stationary distribution is achieved, the ensemble average of 
Eq. (j4"0l) can be approximated by a time average along the CL trajectory. 

In the numerical application of CL dynamics to affect a field-theoretic simulation, physical 
space and time are both discretized. For bulk simulations with periodic boundary conditions 
imposed on the fields, we have found that the spatial discretization is most conveniently 
accomplished by spectral collocation using a plane wave basis and a uniform computational 
grid of M sites |42j. Fast Fourier transforms (FFTs) can then be used to efficiently switch 
between real space (i.e. an M-vector of field values on the grid sites) and Fourier space (i.e. 
the first M Fourier coefficients) representations of the fields. 

Upon spectral collocation, the continuum CL Eqs. (H2|) and (|43l) are transformed into 
a set of 4M stochastic differential equations that can be integrated forward in time from 
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an initial field configuration by standard algorithms 43j, |44j. The simplest algorithm is the 
explicit Euler-Maruyama time integration scheme, 



W {t+At) = W {t) _ Xw At(Ax) 3 + ip®^ + Bj$ 



(44) 



and 



# (*+A*) = # (t) _ A0At(Ax) 3 / + ipc («) j + J ^(*) ( (45) 

where Ax is the grid spacing used for the spacial discretization, and At is the time spacing 
used for the temporal discretization. Here, W, p n , and p c are M-vectors of complex 
variables, which represent the collocated values of the fields W(r), $(r),p n (r), and p c (r) on 
the computational grid. The superscripts index the discrete time at which the M-vectors 
are evaluated. R^ and R^ are M-vectors of real Gaussian random variables, which are 
obtained by spatial collocation of the continuous functions R w (r,t) and R^ir^t) defined by 

/t+At 
drr] w (r,T) (46) 

and 

/t+At 
drri^r). (47) 

The grid-collocated fields R w and R^ are Gaussian white noise functions with vanishing 
averages, (Rw ) = (R^) = 0, second moments given by 

(R^R^) = 2X w At5 tjt/ l (48) 

(« } ) = 2A^A^1 (49) 

and vanishing cross- correlations. The rank M unit tensor is denoted by 1. 

Stochastic time integration algorithms for the CL equations lead to time discretization 
errors in computed expectation values that scale as a power of the time step At. The Euler- 
Maruyama (EM) scheme summarized by Eqs. (jUJ) and (|45|) is of weak order one, which 
implies that the errors in computed averages vanish as At for small At. This low order 
accuracy and the poor stability characteristics of the Euler-Maruyama algorithm make it 
unsuitable for the large-scale 3d simulations reported here. We have instead adopted a 



semi-implicit, weak second order algorithm developed by Lennon and coworkers 45 



was itself inspired by operator splitting methods devised by Petersen and Ottinger 



, which 

3- 



44 



Beyond the second order accuracy, the Lennon algorithm utilizes analytic information about 
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the linearized force in a semi-implicit update scheme to improve stability. The algorithm 
can be cast in a predictor-corrector form with the field updates conducted in Fourier space 
(discrete Fourier transforms of the collocated fields are denoted by carets). The predictor 
steps are explicit EM updates of Eq. (|44|) and Eq. (|45|) in Fourier space: 

\y(*) = W-W - A w + ip%\ + Jig) (50) 

and 

4W = 6W-A^^+i^+A«, (51) 

where the parameters A TO and A<j) are defined as A w = X w At(Ax) 3 and A^ = A^At(Ax) 3 . 
The corrector steps are 

(*) I 6C) 



^ (f+At) = ^ w + [1 + A wPo Ng D (k 2 R 2 g )] - AjpJV + R 

2 + A w , y ww ( K k) 



(52) 



and 

^ t+At) = jj^ + [1 + A^oA^ D (fc 2 ^)] $(*) - A^ip W + gf 

2+ ' 



(53) 

where and Pc*-*- 1 are the segment number density and the charge density operators based 
on the predicted fields of 4> ( *) and . 

This improved stochastic integration algorithm was recently applied to FTS-CL simula- 
tions of block copolymer melts and has been extensively tested in that context by Lennon 



and coworkers 45|. In the present case of polyelectrolyte solutions, we have found that the 
Lennon algorithm permits the use of a time step that is an order of magnitude larger than 
that mandated by stability for the EM algorithm. This translates to a ten-fold reduction in 
computational time. 

In the following section, the FTS results are discussed in the context of the dimensionless 
parameters introduced in Eq. fl36|) : B, C and E. These parameters appear naturally in the 
field theory if all lengths are scaled by the radius of gyration R g , the W field is rescaled 
according to W = NW, and the $ field is rescaled as $ = zN§. With these scalings, 
the CL equations fl50l) - fl53l) for continuous polymer chains depend on the three intensive 
model parameters (B, C, E), the dimensionless simulation cell size L = L/R g , and on three 
dimensionless parameters that relate to the spatial and temporal resolution of the numerical 
algorithm: 

Ax = Ax/R g (54) 
17 



(At) w = X w N 2 At 



(55) 



(At)^ = \ 4> N 2 z 2 At (56) 

For discrete polymer chains, the number of segments iV appears as an additional in- 
dependent parameter in the update equations. However, we shall see that its influence is 
diminished as N is increased to large values approaching the continuous polymer limit. 

At fixed spatial resolution with the Lennon algorithm, we expect second-order conver- 
gence in average properties as the parameters (At) w and (Ai)^ are reduced. All simulation 
data reported below were obtained by setting (At) w and (At)^ to values such that the time 
integration error was negligible compared to the statistical sampling error. In contrast, due 



to ultraviolet divergences in the continuum field theory [28|, |35j , we do not expect conver- 
gence of certain average properties (such as absolute chemical potentials) as the spatial grid 
spacing Ax is refined. However, the structure factors and phase boundaries studied here 
were found to be devoid of ultraviolet divergences, so no special regularization procedures 
were required to isolate the singularities. 



III. RESULTS AND DISCUSSION 



In this section, we summarize and discuss FTS results for the field theory model formu- 
lated in the previous section in the parameter space of B, C, and E. While these three 
parameters completely determine the intensive thermodynamic properties of a system com- 
prised of continuous Gaussian chains, as discussed above, there is an additional independent 
parameter for the discrete Gaussian chains employed in the simulations - the degree of 
polymerization N. We begin with homogeneous solutions of neutral polymers correspond- 
ing to E = 0. The effect of the strength of the excluded volume interaction (B) and the 
degree of polymerization (N) on segment density correlations are studied using FTS, and 
compared with the RPA structure factor. While studying the effect of N, the parameters C 
and B are maintained at constant values. Next, we discuss the effect of the parameter E on 
density correlations and charge correlations in oppositely charged polyelectrolyte solutions 
with fixed B, C, and N. Finally, we construct the phase diagram of our model system in a 
restricted region of the parameter space. Although a more comprehensive study of the phase 
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diagram is worthy to pursue, the high dimensionality of the parameter space and the sig- 
nificant computational requirements both limit the scope of our investigation. Nonetheless, 
our results shall serve to highlight the power and capability of FTS in addressing difficult 
problems in polyelectrolyte complexation. 

The cubic simulation box applied in our simulations has a volume V of (4i? 9 ) 3 (L = 4) 
and is subject to periodic boundary conditions. Occasionally, a larger simulation box with 
V of {8R g ) 3 was used to confirm that finite size effects were not influencing structure or 
thermodynamics. In discretizing the simulation volume, Ax is desired to be smaller than 
the physical length scales of interest (i.e. the correlation lengths for the two fields), but an 
excessively small Ax adds significantly to the computational cost. In most of our simulations, 
Ax was chosen to be R g /S, so the volume of (4i? 9 ) 3 corresponds to a lattice with M = 32 3 
sites. Unlike Ax, which was fixed, At was varied over a broad range to achieve a consistent 
accuracy of time integration depending on the strength of field fluctuations around the 
mean field: (At) Wt $ G [0.001,0.1]. The correlation time is also highly variable throughout 
parameter space. Typically, 10 4 statistically independent field configurations were sampled 
to calculate averages, which corresponds to a number of time steps in the range of 5 x 10 4 
to 5 x 10 5 . The calculation of each structure-factor curve reported here takes from 4 days 
to 7 weeks of CPU time on a single AMD Opteron 248 processor (2.2GHz). However, the 
algorithm scales nearly linearly with domain decomposition across multiple processors, so 
the simulations reported here are ideally suited for a parallel computing environment. 



A. The segment density correlation in neutral polymer solutions 

By setting the charge on the polymer segments to zero (z — E — 0), our model system 
becomes a solution of electrically neutral polymers in a good solvent. This is the Edwards' 
model [l| (Model A of Ref. 28]). When the w field fluctuations are weak, by application 
of Eqs. (1321) and fl36l) . the segment structure factor of Eq. fl38|) can be approximated by the 
dimensionless form, 

^'" h^ < 57 > 

which indicates that the structure factor of segment density is completely determined by 
one parameter, the product of B and C. In deriving this analytical formula, the chain 
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FIG. 1: RPA structure factor for the segment density when BC = 12. Segment structure factors for 
several values of N, calculated from Eq. (I58p . are compared with the result based on the continuous 
Gaussian chain (CGC), cf. Eq. (pp 



molecules were assumed to be continuous Gaussian chains (CGCs), thus the Debye function 
gi){k 2 Rg 2 ) depends only on the unperturbed radius of gyration R g . However, our simulations 
were conducted using discrete Gaussian chains (DGCs) with N segments on each chain, 
for which the corresponding Debye functions depend on N as well as R g . Therefore, the 
corresponding RPA structure factor for DGCs can be written as 

BCg D (k 2 R g 2 ,N) 



u S nn (k,N) 



BCgn(k 2 R g , N) 



where 



38) 



g D (k 2 R g 2 ,N) 



1 



N 



1? ^ 



cxp 



■(kRc 



>2 |» ~ J\ 

' N 



(58) 



(59) 



RPA segment structure factors for DGCs of various are compared with that of the CGC 
in FIG. [TJ At fixed BC, the low k behavior of the structure factor, related to the isothermal 
osmotic compressibility, is independent of the level of discretization of the constituent chain 
molecules. It is only in the high k regime (compared with 27t/R g ) where the discrete nature 
of the DGC model manifests itself by saturating the 1/k 2 decay of the structure factor. 

As expected, the RPA structure factor proves to be a valid approximation to S nn {k) as 
long as the w field fluctuations are weak and approximately Gaussian. This has been verified 
by comparing the RPA structure factor to results obtained from FTS, which incorporates the 
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(a)Dense solution: B — 1 and C = 12 



N= 128 




(b)Moderately-dense solution: _B = 12 and C = 1 

FIG. 2: The segment structure factor in solutions of neutral discrete Gaussian chains: Symbols 
are FTS results calculated from Eq. (|23|) . Lines are RPA structure factors from FIG. [TJ A three 
dimensional lattice of M = 32 3 sites was employed with periodic boundary conditions. The cell 
volume was V = (4i? 9 ) 3 . uq was varied as N was changed to keep B fixed at either 1 or 12. 

full field fluctuation spectrum. In FIG. [2 structure factors derived from FTS for two different 
solution conditions are compared with the RPA structure factor. FIG. [2(a) corresponds to 
the case of a dense solution (C > B) with B = 1 and C = 12, while FIG. [2(b) describes 
a moderately-dense solution with B = 12 and C — 1. However, both share the same RPA 
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predictions of FIG. [I] because the RPA segment structure factor only depends on the product 
BC which is fixed at 12. As FIG. [2(a) shows, there is very good quantitative agreement 
between RPA structure factors and FTS-derived structure factors when the solution is dense 
(C I B = 12). However, in the moderately-dense regime, as exemplified by FIG. [3(b), the 
RPA breaks down, especially at low k. Due to strong excluded volume correlations in this 
regime, the system evidently has a larger osmotic compressibility than is predicted by the 
RPA. 

Additionally, the FTS results at low k show that the strength of the w field fluctuations 
actually depends on N. As implied in Eq. ( |35i) . the w field fluctuations tend to increase 
the osmotic compressibility of the system, and the FTS results indicate that this tendency 
gets stronger for smaller N (see insets of low k regime in FIG. [2]). While a discrete chain 
with iV > 256 is sufficient to model iV-independent thermodynamic properties in a dense 
solution with C/B = 12, a substantially larger iV (> 1024) is required for iV-independent 
thermodynamics in a moderately-dense solution of B/C = 12. 

At larger k and for large N, we can apply an asymptotic expression for the Debye function, 
cjD{k 2 Rg 2 ) ~ 2/k 2 R g 2 , which is highly accurate for kR g > 2tt. Using this approximation, 
the RPA structure factor for the segment density (Eq. (jSTj) ) can be rearranged as 

Q X m ^l + iuk 2 , (60) 
u S nn (k) 

where £ M is the correlation length for segment density. We use Eq. f[60l) to define £ u even 

1 /2 

outside of the RPA, but note that the RPA predicts that £ u = where £ E = R g / (2BC) 
is the Edwards correlation length. In FIG. [31 segment correlation lengths are estimated 
from FTS structure factors, which quantitatively agree with a prior FTS study that used an 



independent method of estimating £ M 



34( | . Not surprisingly, the Edwards correlation length 



is a good approximation to the segment correlation length extracted from FTS when the 
solution is dense. The segment correlation length in a moderately-dense solution, however, 
is underestimated by the Edwards correlation length. 



B. Correlations in symmetric polycation-polyanion mixtures 

In this subsection, density correlations in homogeneous polyelectrolyte solutions are stud- 
ied with FTS. Unlike the results presented so far for uncharged polymer solutions, charge 
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FIG. 3: Inverse structure factor plot used to extract the correlation length for segment density, 
£ u . According to Eq. (|60f) . £ u can be estimated from the ratio of the slope to intercept in the 
intermediate k region of this plot, 1/R g < k/2ir < l/£ u . We find that £ U /R g = 0.201 ± 0.0014 in 
the case of B = 1 and C = 12, and £ U /R g = 0.225 ± 0.0016 for B = 12 and C = 1. The Edwards 
correlation length is ^E/Rg — 0.204 for BC = 12. The simulations were conducted on a M = 32 3 
lattice with periodic boundary conditions, a system volume of V = (4R g ) 5 , and a chain length of 
N = 1024. The linear fit was applied in the regime of {kRg) z £ [30,200]. 

correlations, as well as segment correlations, are of interest in polyelectrolytes. Specifically, 
we are interested in the interplay between charge correlations and segment correlations and 
how they relate to the phenomenon of complex coacervation. In all the FTS results reported 
in this subsection, the solution was dense (C = 12, B — 1) and discrete Gaussian chains of 
iV = 1024 were used. 

The RPA structure factor of Eq. (1571) for the segment density is evidently independent of 
the charge content in the system; this is a result of decoupling of w and </> fluctuations at the 
RPA level. It would seem that without a higher-order analysis of fluctuations, it is impossible 
to guess even the qualitative effect of charge correlations on segment correlations. However, 
the one loop result for the osmotic pressure does provide some insight. From Eq. (15BT) . the 
charge correlations are seen to reduce the osmotic pressure of a homogeneous solution, the 
same effect caused by the segment correlations. Thus, by adding equal and opposite charges 
to the chains of a neutral polymer solution, we can expect to obtain a solution with increased 
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FIG. 4: The segment structure factor in solutions of symmetric polyelectrolytes. The symbols are 
FTS results and the line is calculated from the RPA result, Eq. (I58p . The following parameters 
were used: B = 1, C = 12, N = 1024, V = 4 x 4 x AR g 3 , and M = 32 x 32 x 32. 

osmotic compressibility. This increase in compressibility will in turn be manifest in the low 
k behavior of S nn (k). 

In FIG.IU segment structure factors for symmetric polyelectrolyte solutions were obtained 
from FTS. Because the solution is dense, the RPA structure factor is observed to be a good 
approximation when the polymers are neutral (E = 0). However, as the charge density is 
increased from zero, the FTS results show strong deviations from the RPA structure factor 
which is independent of E. The osmotic compressibility increases monotonically with E and 
eventually diverges at even higher E, indicative of a macrophase separation (in this case 
"complex coacervation" ) . 

Another structure factor of interest in the charged system relates to the charge density 
correlation. When the <fi field fluctuations are weak, the RPA formula for the charge density 
structure factor, Eq. (J55|) . can be combined with Eqs. (TH3T) and (I5rj|) to obtain the following 
approximation: 

**lB q (h) „ C Egp{k 2 R g 2 ) 
k 2 cc[ >~ k 2 R 2 g + CEg D (k 2 R g 2 Y 1 J 

This RPA formula also applies to solutions of discrete Gaussian chains when the modified 

Debye function of Eq. (|59|) is substituted for Eq. (!30l) . We see from Eq. (I6T1) that the RPA 

structure factor for the charge density is completely dictated by the combination parameter 

24 




FIG. 5: The charge structure factor in solutions of oppositely charged polyelectrolytes. Symbols are 
FTS results, while the lines are RPA structure factors calculated from Eq. (I6ip . In the simulations, 
we vary the parameter E by changing z while keeping l B constant. The parameters correspond to: 
l B = Rg/8, B = 1, C = 12, N = 1024, V = (4R g ) 3 , and M = 32 3 . 

CE, and is independent of B. 

By again utilizing the asymptotic expression for the Debye function, g B (k 2 R g 2 ) « 
2/k 2 R g 2 , valid for N ^> 1 and kR g > 2n, the inversion of Eq. (!61~!) provides 

7 , ^ 7T\ « 1 + e c ^ 4 , (62) 

which defines a new charge correlation (or screening) length £ c . Explicit use of the RPA 
formula leads to the result £ c = £p#, where ^pe = R g /(2EC) 1 ^ . 

Eqn (|62l) also predicts that S cc (k) is maximal at = l/£ c . Note that the correlation 
length Cpe is proportional to the —1/4 power of the segment density, and hence is quali- 
tatively different from the Debye-Hiickel length for small ions, which is proportional to the 
— 1/2 power of the ion density. Thus, the attachment of charges to polymer chains cre- 
ates a coupling between chain conformational statistics and charge density that drastically 
changes the electrostatic correlation properties of polyelectrolyte solutions compared with a 
conventional small-ion electrolyte. 

In FIG. El FTS results for the charge density structure factor are compared with RPA 
predictions based on Eq. (16TI) . The dimensionless object 4iil B Rg 2 S cc , instead of Anl B S cc /k 2 , 
is plotted to clearly show the location of the maximum (~ l/£ c ) an d the screening of charge 
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density at low k: S cc (k) ~ k 2 . Although the RPA remains a good approximation at low 
charge content (E = 900, 3600), the deviation between the RPA and FTS becomes no- 
ticeable at larger E and larger k. As E increases, the maxima of both the RPA and FTS 
results shift to higher k, indicating a shorter screening length consistent with the RPA scal- 
ing £ c ~ (EC)~ l l 4 . At E = 19600, however, the RPA structure factor places the maximum 
at a somewhat higher value of k than FTS. In other words, the RPA underestimates the 
charge correlation length for strongly charged chains. It is also evident that the RPA slightly 
overestimates the amplitude of charge correlations at large E. Later, it will be shown that 
this qualitative observation is consistent with differences observed between phase boundaries 
deduced from FTS studies and the one-loop analysis. 



C. Complex coacervation: phase diagram 

As implied by the diverging osmotic compressibility in FIG. HJ a macrophase separation 
(complex coacervation) is possible in our field theory model. It was already anticipated in 
the analytic one-loop correction to the osmotic pressure, Eq. (1351) . that a phase separation 
may result from the competition between the positive second virial term from the excluded 
volume interaction and the negative contribution from charge correlations. The spinodal 
(single-phase stability limit) is determined by the condition d/3p/dC = from the one loop 
osmotic pressure expression of Eq. ( 1351) . With an additional approximation regarding the 
dilute phase, the one-loop theory can also predict the binodal (the coexistence curve of dilute 
and coacervate phases). Assuming a supernatant dilute phase free of polyelectrolyte, we can 
write the binodal equation in the form of (3p = 0. 

In the context of computer simulations, a rigorous construction of a phase diagram re- 
quires a computational technique for accessing the free energy of the system. In the case of 
particle-based simulations, a variety of free energy estimation methods are available, includ- 



ing thermodynamic integration 
and histogram techniques [47 . 



;echniques, Gibbs ensemble and particle insertion methods, 



481 ]. Such methods, however, are only now emerging for 
field-based simulations, so here we make a crude estimate of the boundary of the two-phase 
region in our polyelectrolyte model by monitoring the hysteresis of an order parameter while 
varying an intensive parameter of the system. The order parameter chosen to characterize 
complex coacervation is the density difference Ap between the dense coacervate phase and 
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FIG. 6: Phase diagram for the three-dimensional symmetric polyelectrolyte mixture expressed 
in the coordinates of reduced polymer concentration C and reduced excluded volume B at fixed 
reduced Bjerrum length E = 14400. Solid and dashed lines are the analytical one-loop binodal 
and spinodal, respectively. The one-phase region is above and to the right of the lines, and the 
two-phase region is below and to the left. Symbols are from the examination of the hysteresis 
in Ap using FTS; details of the procedure are explained in the text. FTS data obtained from 
supercooling simulations (lowering B) are denoted by •, while superheating data are denoted by 
■. The dotted line is a power-law fit. The error in the numerical data is comparable to the size 
of the symbols. Numerical simulations were conducted in a cubic cell of size (4R 9 ) 3 with periodic 
boundary conditions and with iV = 384. 

the dilute phase. Hysteresis is examined while varying the solvent quality B. For example, 
consider a homogeneous single phase solution (Ap = 0) at certain values of B, C, and E. 
While slowly "cooling" that system by gradually decreasing B at fixed C and E, a sudden 
jump in Ap from zero to a finite value occurs as the system separates into two phases of 
different densities. The high density phase can be identified as the complex coacervate. On 
the other hand, upon "heating" the system from the two-phase region by gradually increas- 
ing B at fixed C and E, a sudden drop of Ap from a finite value to zero is observed as 
the system exceeds the limit of superheating and remixes into one homogeneous phase. For 
most phase transitions, the phase coexistence curve is closer to the superheating curve than 
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to the supercooling curve. Thus, while the superheating and supercooling curves should 
bracket the transition, we expect that the superheating curve will lie closer to the binodal 
boundary. 

An example of a phase diagram (in 3D) constructed in such a way is provided in Fig. [HI 
Spinodals and binodals are surfaces in the three-parameter space of the reduced variables 
C, B, and E. The figure represents a cross-section of this three-dimensional space by a 
plane E = 14400; hence, the diagram involves only the C and B variables. The one-phase 
region (disordered homogeneous phase) is above and to the right of the lines, and the two- 
phase region is below and to the left. The tie lines in the two-phase region are horizontal 
(constant B) and connect nearly pure solvent (C ~ 0) with a coacervate phase at the binodal 
concentration. Remarkably, the analytical and numerical results nearly coincide in the high 
concentration region of the figure, despite the limitations of each method. The numerical 
results are subject to finite cell size and chain discretization limitations (N = 384), while the 
analytical predictions neglect two-loop and higher order terms in the fluctuation analysis. 
Nonetheless, our numerical supercooling result practically follows the analytical spinodal, 
and the analytical binodal yields similar exponent (—1.31) as is obtained from a power-law 
fit to the numerical superheating points (—1.38). The over-estimate of the size of the two 
phase region by the one-loop theory is consistent with the observation in FIG. [5] that the 
RPA structure factor over-estimates the strength of charge correlations at large E and hence 
expands the two-phase region. The overall semi- quantitative agreement between theory and 
FTS, however, indicates that both approaches have utility for this class of problems. 



IV. CONCLUSIONS AND PERSPECTIVE 



In this paper, we reported on the application of the emerging field-theoretic computer 
simulation (FTS) technique to a simple model of polyelectrolyte complexation. Specifically, 
we built and numerically simulated a field theory model of a salt-free solution containing a 
symmetric mixture of flexible polyanions and polycations in a good solvent. This particular 
system constitutes a minimal model for the phenomenon of complex coacervation, a type of 
liquid-liquid phase separation in which a nearly pure solvent phase coexists with second fluid 
phase (the "coacervate") that contains the majority of the polyelectrolytes. Theoretically, 
the symmetric coacervate model is interesting because the workhorse tool of polymer physics, 
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self-consistent field theory (SCFT), fails to describe the electrostatic effects responsible for 
coacervation. 

Previous analytical work on closely related models of polyelectrolyte mixtures has shown 
that complex coacervation can be predicted based on calculations that assume weak, 
Gaussian field fluctuations, i.e. calculations at the one-loop level of fluctuation expan- 
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261 ] . However, the reliability of these predictions have remained 



unclear because analytical techniques for treating more realistic situations of strong charge 
and excluded volume correlations are lacking. The FTS results of our paper are significant 
because they can provide numerical data to assess the validity of the RPA, both in terms 
of its predictions for charge and segment density correlations and for the location of the 
two-phase envelope. Overall, we have found that the RPA is remarkably robust in its pre- 
dictions, except at very high charge densities (large values of the parameter E) where it 
overestimates the strength of charge correlations and the size of the two-phase region. 

Perhaps more significantly, our results have validated the emerging field theoretic polymer 
simulation technique as a powerful new tool for examining the structure and thermodynamics 
of polyelectrolyte systems. FTS can be applied even in situations where the RPA is inap- 
plicable or technically very difficult, such as cases of mesophases formed by weakly charged 
polymers with hydrophobic backbones or block co-polyelectrolytes [37]. Complexation of 
charged polymers that also contain neutral blocks or grafts (which can be either hydropho- 
bic or hydrophilic) can also produce inhomogeneous "structured coacervate" phases that 

nnnn 

are not amenable to study by current theoretical methods [15|, H6|, 113, L18|] . We expect that 
FTS will prove to be a valuable tool for exploring these and related types of polyelectrolyte 
systems. 
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FIG. 1: RPA structure factor for the segment density when BC = 12. Segment structure 
factors for several values of N, calculated from Eq. fl58l) . are compared with the result based 
on the continuous Gaussian chain (CGC), cf. Eq. (1BT1) . 

FIG. 2: The segment structure factor in solutions of neutral discrete Gaussian chains: 
Symbols are FTS results calculated from Eq. (1231 . Lines are RPA structure factors from 
FIG. [TJ A three dimensional lattice of M = 32 3 sites was employed with periodic boundary 
conditions. The cell volume was V = (4i? ff ) 3 . u was varied as N was changed to keep B 
fixed at either 1 or 12. 

FIG. 3: Inverse structure factor plot used to extract the correlation length for segment 
density, £ u . According to Eq. (1601) . £ u can be estimated from the ratio of the slope to 
intercept in the intermediate k region of this plot, 1/R g < k/2ii < l/£ u . We find that 
i u /R g = 0.201 ± 0.0014 in the case of B = 1 and C = 12, and £ U /R g = 0.225 ± 0.0016 
for B = 12 and C = I. The Edwards correlation length is £ E /Rg - 0-204 for BC = 12. 
The simulations were conducted on a M = 32 3 lattice with periodic boundary conditions, a 
system volume of V — (4R g ) 3 , and a chain length of N = 1024. The linear fit was applied 
in the regime of {kRgf G [30,200]. 

FIG. 4: The segment structure factor in solutions of symmetric polyelectrolytes. The 
symbols are FTS results and the line is calculated from the RPA result, Eq. ( 1581) . The follow- 
ing parameters were used: B = 1, C = 12, N = 1024, V = 4 x 4 x 4R g 3 , and M = 32 x 32 x 32. 

FIG. 5: The charge structure factor in solutions of oppositely charged polyelectrolytes. 
Symbols are FTS results, while the lines are RPA structure factors calculated from Eq. (I6T!) . 
In the simulations, we vary the parameter E by changing z while keeping Ib constant. The 
parameters correspond to: l B = R g /8, B = 1, C = 12, iV = 1024, V = (4i? 9 ) 3 , and M = 32 3 . 

FIG. 6: Phase diagram for the three-dimensional symmetric polyelectrolyte mixture 
expressed in the coordinates of reduced polymer concentration C and reduced excluded 
volume B at fixed reduced Bjerrum length E = 14400. Solid and dashed lines are the 
analytical one-loop binodal and spinodal, respectively. The one-phase region is above and 
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to the right of the lines, and the two-phase region is below and to the left. Symbols are from 
the examination of the hysteresis in Ap using FTS; details of the procedure are explained 
in the text. FTS data obtained from supercooling simulations (lowering B) are denoted by 
• , while superheating data are denoted by ■. The dotted line is a power-law fit. The error 
in the numerical data is comparable to the size of the symbols. Numerical simulations were 
conducted in a cubic cell of size (4i? ff ) 3 with periodic boundary conditions and with N = 384. 
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